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Abstract. I propose an iterative expectation maximization algorithm for 
reconstructing the density matrix of an optical ensemble from a set of balanced 
homodyne measurements. The algorithm applies directly to the acquired data, 
bypassing the intermediate step of calculating marginal distributions. The 
advantages of the new method are made manifest by comparing it with the 
traditional inverse Radon transformation technique. 
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Quantum tomography is a technique of characterizing a state of a quantum system 
by subjecting it to a large number of quantum measurements, each time preparing 
the system anew. By varying the configuration of the measurement apparatus, one 
acquires the quantum statistics associated with different bases from which complete 
information about the state of the system can be extracted. 

The ensemble's density matrix can be evaluated from the experimental statistical 
data by a number of techniques. In this paper we are dealing with one such technique, 
the Maximum Likelihood (MaxLik) estimation. Assuming a particular density matrix 
p, one can evaluate the likelihood (probability) of acquiring a particular set of 
measurement results. The ansatz of the MaxLik method is to find, among the variety 
of all possible density matrices, the one which maximizes the probability of obtaining 
the given experimental data set. To date, this method has been applied to various 
quantum and classical problems from quantum phase estimation ^ to reconstruction 
of entangled optical states [31 

In the present work we discuss applications of likelihood maximization to quantum 
homodyne tomography - a method of characterizing a quantum state of an optical mode 
by means of multiple phase-sensitive measurements of its electric field's quantum 
noise. Since its proposal |31 Eind first experimental implementation |S] in the early 
1990s, quantum homodyne tomography has become a robust and versatile tool of 
quantum optics and has been applied in many different experimental settings. As 
any statistical method, it is compatible with the likelihood maximization approach. 
Yet so far experimentalists have used a clumsier and less accurate reconstruction 
algorithms based on the inverse Radon transformation. The goal of this paper is to 
provide an explicit numerical recipe for the MaxLik reconstruction of a quantum state 
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from homodyne data and to demonstrate its advantages by applying it to bona fide 
experimental results. 

The applications of MaxLik estimation to homodyne tomography have been 
investigated by Banaszek, who has reconstructed the photon-number distribution (the 
diagonal density matrix elements which correspond to a phase-randomized optical 
ensemble) from a Monte-Carlo simulated data set jHlIZ]. In a subsequent publication 
|H] , Banaszek et al. discussed the MaxLik estimation of the complete density matrix, 
but no explicit reconstruction algorithm has been presented. 

The iterative scheme outlined here is based on that elaborated by Hradil et al. 
for discrete- variable states [SJ IIOL [TT] . Here we give its brief overview. Consider a 
large set of von Neumann measurements, each one projecting the state of the system 
onto an eigenstate of a measurement apparatus \yj). The set of all possible outcomes 
can be associated with either one or several measurement bases. Let fj be 
the frequency of occurrences for each outcome. Then, with the system being in the 
quantum state p, the likelihood of a particular data set {fj} is 

C{p)^X{w^. (1) 

where pr^ = {yj\p\yj) — Trpjp] is the probability of the outcome and lij ~ 
\yj) {yj\ denotes the projection operator. 

To find the ensemble p which maximizes the likelihood Q]), one introduces the 
operator 

i?(p) = E— (2) 

and notices that for the ensemble po that is most likely to produce the observed data 
set, fj oc pr^. Furthermore, since the lij oc 1, we find R{po) oc 1 and thus 

R{po)po PaR{M « Pa (3) 

as well as 

R{po)poR{po) oc po. (4) 

The last relation forms the basis for the iterative algorithm. We choose some 
initial denstity matrix as, e.g., p^"-* = A/'[i], and apply repetitive iterations 

= N , (5) 

where Af denotes normalization to a unitary trace. Each step will monotonically 
increase the likelihood associated with the current density matrix estimate while 
the latter will asymptotically approach the maximum-likelihood ensemble po* . This 
iterative scheme can be viewed as a special case of the classical expectation- 
maximization algorithm |12| . 

Now we turn to the main subject of the paper and consider a homodyne 
tomography experiment performed on an optical mode prepared in some quantum 
state p. In an experimental run one measures the field quadrature at various 
phases of the local oscillator. Each measurement is associated with the observable 
Xg = X cos9 + PsinO, where X and P are the canonical position and momentum 
operators and 6 is the local oscillator phase. 

* We base our iteration scheme on Eq. jlj rather than ^ in order to ensure the positivity of the 
density matrix at each step. 
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For a given phase 0, the probabiUty to detect a particular quadrature value x is 
proportional to 

pvg{x)^TT[Ili9,x)pi (6) 

where 11(6', x) = \6,x){6,x\ is the projector onto this quadrature eigenstate. In the 
Fock (photon number state) basis, the projection operator is expressed as 

TimniO, x) = (m| n(0, x) \n) = (m|0, x) (0, x\n) , (7) 

where the overlap between the number and quadrature eigenstates is given by the well 
known stationary solution of the Schrodinger equation for a particle in a harmonic 
potential: 

{n\0,x)^e-<^ f^Y'^h^ exp(-x2), (8) 

with Hn denoting the Hermite polynomials*. 

Because a homodyne measurement generates a continuous value, one cannot apply 
the iterative scheme © directly to the experimental data. One way to deal with this 
difficulty is to discretize the data by binning it up according to 6 and x and counting 
the number of events fe,x belonging to each bin. In this way, a number of histograms, 
which represent the marginal distributions of the desired ensemble's Wigner function, 
can be constructed. They can then be used to implement the above reconstruction 
procedure. 

However, discretization of continuous experimental data will inevitably lead to a 
loss of precision. To lower this loss, one needs to reduce the size of a single bin and 
increase the number of bins. In the limiting case of infinitely small bins, fe^x takes on 
the values of either or 1, so the likelihood of a data set {{6i,Xi)} is given by 

In/: = ^lnprg^(x,), (9) 

i 

and the iteration operator |(5J) becomes 



£,,.x n(g„x,) 



where i = 1 . . . N enumerates individual measurements. The iterative scheme ((SJ can 
now be applied to find the ensemble which maximizes the likelihood ©. 

In practice, the iteration algorithm is executed with the density matrix in the 
photon number (Fock) representation. Since the Hilbert space of optical states is of 
infinite dimension, the implementation of the algorithm requires its truncation so the 
Fock terms above a certain threshold are excluded from the analysis. This assumption 
conforms to many practical experimental situations in which the intensities of fields 
involved are a priori limited. 

It is instructive to compare the maximum likelihood quantum state estimation 
with the traditional methods employed in homodyne tomography: reconstruction of 
the Wigner function by means of the inverse Radon transformation )13| and evaluation 

* Normalization [x,p\ = i/2 is used. The additional phase factor e'™^ originates from the properties 
of the phase-space rotation operator 1131 U{9) = e~'^". From U'^ (0)aty{9) = ae^'* we find for 
the quadrature operator IJ^ {9)XU{9) = Xg and for its eigenstate \9,x) = \0,x). From the first 
and last relations above, we obtain {m\9,x) = e*^"* {m\0,x). The quantity {m\0,x) is the energy 
eigenwavefunction of a harmonic oscillator. 




Figure 1. Estimation of an optical ensemble from a set of 14152 experimental 
homodyne measurements 1161 by means of the inverse Radon transformation (a) 
and the likelihood maximization algorithm (b). The Wigner function and the 
diagonal elements of the reconstructed density matrix are shown. The inverse 
Radon transformation in (a) was performed by means of the filtered back- 
projection algorithm with the cutoff frequency of 6.3. The statistical uncertainties 
in (b) were determined by means of a Monte-Carlo simulation (see text). 



of the density matrix using quantum state sampling 1141 115j . Fig. ^ shows apphcation 
of these two techniques to the experimental data from Ref. |TC| . The data set consists 
of 14152 quadrature samples of an ensemble approximating a coherent superposition 
of the single-photon and vacuum states. 

The reconstruction shown in the figure reveals the advantages of the MaxLik 
technique in comparison with the standard algorithm. First, the finite amount and 
a discrete character of the data available leads necessarily to statistical noise which 
prevents one from extracting complete information about a quantum state of infinite 
dimension. To deal with this issue, both techniques apply certain assumptions on 
the ensemble to be reconstructed. While the MaxLik algorithm truncates the Fock 
space, the filtered back-projection imposes low pass filtering onto the Fourier image 
of the Wigner function*, i.e. assumes the ensemble to possess a certain amount of 
"classicality" ^Tj. The latter assumption is dictated by mathematical convenience and 
is much less physically founded than the former. The ripples visible in the Wigner 
function reconstruction in Fig. ^a) are a direct consequence of statistical noise and 
are associated with unphysical high number terms in the density matrix. Such ripples 
are typical of the inverse Radon transformation ^1 . 

* The pattern function reconstruction of the density matrix is free from this drawback as it does not 
involve spectral filtering and relies on truncating the Fock space instead. 
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Second, the back-projection algorithm does not impose any a priori restrictions on 
the reconstructed ensemble. This may lead to unphysical features in the latter, such as 
negative diagonal elements of the density matrix in Fig. ^a). The MaxLik technique, 
on the other hand, allows to incorporate the positivity and unity-trace constraints 
into the reconstruction procedure, thus always yielding a physically plausible ensemble 

lElIZI 

A third important advantage of the MaxLik technique is the possibility to 
incorporate the detector inefBciences. In a practical experiment, the photodiodes 
in the homodyne detector are not 100% efficient, i.e. they do not transform every 
incident photon into a photoelectron. This leads to a distortion of the quadrature 
noise behavior which needs to be adjusted for in the reconstructed ensemble. 

A common model for a homodyne detector of non-unitary efficiency 77 is a perfect 
detector preceded by a fictitious beam splitter of transmission 77. The reflected mode is 
lost so the incident density matrix undergoes, in transmission through the imaginary 
beam splitter, a so-called generalized Bernoulli transformation |13| : 

00 

(m| pr, \n) = ^ B,n+k,m{Tl)Bn+k,n{'n) {m + k\ Pq \n + k) , (11) 

k=Q 

where po and p,j are the density matrices of the original and transmitted ensembles, 
respectively, and Bn+k.n = \j {^ti^^ ?7"(1 — 77)*-'. Under these circumstances the 
probability © of detecting a quadrature value x becomes 

w'lix)^ {e,x\p^\e,x) (12) 

00 00 

7i-\-k.'n 

(rj) {n\6, x) {9, x\m) {m + k\ po \n + k) , 

m,n— k—0 

SO the projection operator 11(6', x) becomes replaced by a POVM element given by 
En{9,x)^ B,n+k.rn{v)Bn+kA:n){n\0,x) {9,x\m)\n + k) {vi-^k\. (13) 



ni.n 



k 



Aware of the homodyne detector efficiency rj, one runs the iterative algorithm |(SJl and 
reconstructs the original density matrix po (S]. 

Theoretically, it is also possible to correct for the detector inefficiencies by 
applying the inverted Bernoulli transformation after an efficiency-uncorrected density 
matrix has been reconstructed ,20). However, this may give rise to unphysically large 
density matrix elements associated with high photon numbers. With the inefficiency 
correction incorporated, as described above, into the reconstruction procedure, such 
terms do not arise [7]. 

It is interesting to discuss the MaxLik reconsruction of the density matrix in 
comparison with the point-by-point reconstruction of the Wigner function as proposed 
in 1511 . To determine the value of the Wigner function at a specific point {p, q) in the 
phase space, Ref. 2L proposes to apply a phase-dependent shift to the experimental 
data which corresponds to a displacement of the point (p, q) to the phase space origin. 
Then one reconstructs a phase-averaged ensemble according to Refs. U\, and 
calculates the value of the Wigner function at the origin of the displaced phase space, 
which is equal to that at the desired location in the undisplaced phase space. 

This scheme may appear more involved than the one proposed here, as one needs 
to run a separate iteration series for every point in which the Wigner function is to be 
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calculated. However, due to a smaller number of parameters and a simplified iteration 
step, each iteration takes less time and the series converges faster. The choice of a 
particular scheme depends on a specific task and on the chosen truncation threshold 
in the Fock space. It is important to note that the scheme j21j does not impose any a 
priori restrictions onto the ensemble being reconstructed, and therefore the latter is 
not guaranteed to be physically meaningful. 

Finally, we discuss statistical uncertainties of the reconstructed density matrix. 
In generic MaxLik algorithms, they are typically estimated as an inverse of the Fisher 
information matrix [221113 G = d'^ C{t) / dtdV , where t denotes a set of independent 
parameters with respect to which the likelihood is evaluated. Because the density 
matrix elements are not fully independent but bound by the positivity and unity 
trace constraints, one expresses the density matrix as a product p = ft7^/Tr[ftf]. 
Now the constraints for p are satisfied for any random T and one can regard the 
elements of the latter as free parameters in evaluating the Fisher information O . 

A sensible alternative is offered by a clumsy, yet simple and robust technique of 
simulating the quadrature data that would be associated with the estimated ensemble 
Pml if it were the true state. One generates a large number of random sets of 
homodyne data according to Eq. ®, then applies the MaxLik reconstruction scheme 
to each set and obtains a series of density matrices p'^. each of which approximates the 
original matrix pml- The average difference {\pml — p'k\)k evaluates the statistical 
uncertainty associated with the reconstructed density matrix. 
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